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Abstract 

We discuss the efficient implementation of Hamiltonian BVMs (HBVMs), a recently introduced class 
of energy preserving methods for canonical Hamiltonian systems (see [2| and references therein) , via 
their blended formulation. We also discuss the case of separable problems, for which the structure 
of the problem can be exploited to gain efficiency. 
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1. Introduction 

The conservation of energy allows to avoid the numerical drift observed when using standard 
numerical methods for solving canonical Hamiltonian problems, i.e., problems in the form 

y' = JVH(y), J=(° Im 7 o n )> i/(*o) = Vo e M 2m , (1) 

where H(y) is a smooth scalar function and, in general, I r will hereafter denote the identity matrix 
of dimension r (when the lower index will be omitted, the size of the matrix can be deduced from the 
context). In this respect, Hamiltonian Boundary Value Methods (HBVMs) is a recently introduced 
class of methods able to conserve energy when H(y) is a polynomial of arbitrarily high degree. 
Clearly, this implies a practical conservation of energy for any suitably regular Hamiltonian function, 
which will be assumed hereafter. We refer to @, S S S 1 an d references therein, for an overview 
on energy-conserving methods and the derivation of HBVMs. When problem ([1]) is separable, i.e., 
when 

H( y )=H(q,p) = ^p T p-U(q), q,peR m , (2) 
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then JT]) reduces to a special second order equation, 

q" = VU(q), 

and the associated HBVM may be properly formulated in order to take advantage, in terms of 
efficiency, from the above simplification. 

In this paper we investigate the efficient implementation of HBVMs, also in the case of sepa- 
rable problems. In more details, in Section [2] we briefly derive HBVMs. Then, in Section [3] we 
investigate the efficient solution of the generated discrete problem, via the blended implementa- 
tion of the methods, which has already proved to be very effective in other settings (sec, e.g., 
P, |H, 0, [13, EH El , 13 , 14 [ ) . The case of separable problems is then discussed in Section 2J A few 
numerical tests, along with some concluding remarks are then given in Section [5l 



2. Hamiltonian BVMs (HBVMs) 

The derivation of HBVMs will be done according to the approach followed in @, 0] , which further 
simplifies the already simple idea initially used in 

HHHH (see also (IEG3). Let 

us then considcr 

the restriction of problem ([TJ to the interval [to, to + h], with the right-hand side expanded along an 
orthonormal basis {Pj}j>o : 



'(t +rh) = J^P j (T) f P 3 {c)VH{y{to+ch))&c, r € [0, 1] 
j>0 Jo 



(3) 



In particular, we here consider an orthonormal polynomial basis, provided by the shifted and scaled 
Legendre polynomials on the interval [0, 1], even though the arguments can be easily extended to 
more general bases. The basic idea, is now that of looking for an approximate solution belonging to 
the set of polynomials of degree not larger than s. This is achieved by truncating the series at the 
right-hand side in ([3]), thus obtaining the approximate problem 



s-l -i 

a'(t + Th) = jy^P j (T) Pj{c)VH{a{to+ch))&c, rg [0,1], <r(t ) = 

3=0 J ° 



yo- 



(4) 



The approximation to y(to + h) is then given by 

yi =a(t + h). (5) 
The method can be easily seen to be energy-preserving since, considering that J is skew-symmetric, 



H( yi )-H(y ) = h [ \7H(a(t +Th)) T a'(t Q + Th)dT 
Jo 



s-l r i 



3=0 





T 






tIi)) At 


J 


[ Pj-(c)Vfr(c7-(io + cft))dc 
Jo 


= 



Integrating both sides of the first equation in ([4]) yields 

S — 1 pi 

(r{t Q + Th) = y + p j( x )d x Pj{c)JVH(a(t + ch))dc, 

j=o J° J° 



(G) 
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which may be exploited to determine the shape of the unknown polynomial a, provided that a 
technique to handle the rightmost integrals is taken into account: the obvious choice is the use 
of quadrature formulae. If we assume that H(y) is a polynomial of degree v, then the integrals 
appearing in (0| can be exactly computed by a Gaussian formula with k abscissas {cj}, in the event 
that 

k > su/2, (7) 
thus obtaining a discrete problem in the form 

S — 1 « r ■ 



a(t Q + ah) = a l 



, k, 



(8) 



where the \b{\ are the quadrature weights of the formula defined over the abscissae {c{\. For general, 
suitably regular (e.g., analytical) Hamiltonian functions, we can still use formula © in place of 
(jSJ), provided that the integrals in © are approximated to machine prccisiorOJ in the following, 
we will always assume such an accuracy level when a non polynomial function is considered, and 
consequently we will make no distinction between the integrals and the corresponding approximations 
as well as between the two polynomials a obtained by solving either © or ([5]) (see Q for more 
details). 

Method ©-© is called HBVM(fc,s): it was shown 3 that its order is 2s, for all k > s. In 
particular, for k = s it reduces to the well known s-stages Gauss method. 

By introducing the matrices = diag(&i, . . . , bk) and 



T-s-x = 



-Pj_i(x) dx 



P k. X 8 



i = 1 . . . k 

= 1... 8 



V r -1 = (Pf-lte)) l = 1 . 

J = I- 



fkxr 



the HBVM(fc,s) can be recast as a Runge-Kutta method with Butcher tableau 



Cl 








Cfc 






h ••• b k 



(9) 



The next result follows from well-known properties of Legendre polynomials (hereafter e, denotes 
the ith unit vector in W). 



Lemma 1. 



whe 



X„ = 



( 5 

ii 
\ 



s-l 



VsX s = V s 



-6 





-6-1 
o / 



i > 1. 



(10) 



(11) 



1 As we will sec, increasing the order of the quadrature formula, namely the integer k, will not result into an increase 
of the computational cost associated with the implementation of the method. 
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Consequently, the matrix in the Butcher tableau © can be written as 



(12) 



Notice that, since V S X S has s linearly independent columns, the k x k coefficient matrix A has 
rank s: it is then possible to recast the discrete problem in a more convenient form, which clearly 
shows that its (block) size is s, rather than k (see also Q). For this purpose, let us define the (block) 
vectors (see g]) and ©) 

k 

7j = hPj(c e )JVH(a(t + c e h)) : j = 0,...,s-l. (13) 









f 70 \ 






— 






\ k ) 







In view of ((4]) , we see that the vectors jj may be interpreted as the coefficients in the expansion 
of the degree s — 1 polynomial a'(t + rh) along the orthonormal basis {Pj}j=o,..., s -i- 

From ([8]) one obtains 

y = e ® ?/o + /iZa-i ® /zm 7, (14) 
with e=(l,...,l) T eR fc , and then, by virtue of (fT3"]) . one has to solve the equation in the unknown 
7 

F(l) = 7 - CPj-i^ ® J) Vff (e ® y + ® hml) = 0. (15) 

The application of the simplified Newton iteration for solving (|15[) yields 

[I - hV^nis.! ® Go] A £ = -F (V ), Y +1 =Y + A', (16) 

with Go = (j\7 2 H(yo)). By virtue of (fT0|) . and considering that 

Vf^ttPa = (Is 0) £l ,X!+1 , (17) 

(HHJ) reduces to 

[/-/iX s ®Go]A^ = -F( 7 ^), 7 f+1 = 7 ^ + A f , £ = 0,1,..., (18) 
which, as is readily seen, has (block) size s, rather than k. 



3. Blended implementation 

From the arguments in the previous section, one then concludes that the discrete problem, to 
be solved at each integration step when approximating the Hamiltonian problem (TTJ), is given by 
(|T5|) . thus requiring the solution of ([TBI) . We are going to solve such equation by means of a blended 
implementation of the method, according to El EH- Indeed, such implementation of block 

implicit methods has proved to be very effective, leading to the development of the codes BiM 9( 
and BiMD [l3| for stiff ODE IVPs and linearly implicit DAEs (the codes are available at the url fl7|). 
Let us, for sake of simplicity, discard the iteration index. Consequently, we have to solve the linear 
system 

(I - hX s ® Go) A = —F(~f) = rj. (19) 
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Considering that matrix X s (see (jlllO is nonsingular, such equation can be equivalently written as 
p (X; 1 ® hm - hl s ® G ) A = pX- 1 ® J 2m »7 = tji, (20) 
where p is a positive constant. By introducing the (matrix) weight function 

= / s ®£o\ S = (I 2 m-p/iGo)" 1 , (21) 

we then obtain the following problem, which has still the same solution as (fT9"| : 

T(A) = 9 [{I- hX s <8 G ) A - rj] + (I - 9) [p <8 7 2m - W s ® Go) A - m ] = 0. (22) 

One easily realizes that it is obtained as the blending, with weights and {1—6), of the two equivalent 
problems (fTi?)) and (|2()p. respectively. Problem (|22p defines the blended method associated with the 
original one, which we call blended HBVM, in the present case. The free parameter p is chosen in 
order to optimize the convergence properties of the corresponding blended iteration, 



A n+1 = A n - 9T{A n ), 



n > 0, 



(23) 



with an obvious meaning of the lower index. Such iteration only requires (see (1211) ) the factorization 
of the matrix Eo having the same size as that of the continuous problem. According to the linear 
analysis of convergence in (llj . the free parameter p is chosen as 



min{|A| : A e (t{X s )}, 



(24) 



which provides optimal convergence properties (in particular, an L-convergent iteration (ll|V A few 
values of ([Ml) are listed in the table below, for sake of completeness. 



0.2887 0.1967 0.1475 0.1173 



Remark 1. A nonlinear version of i23\) can be readily derived, by taking A n = and updating the 
vectors r) and r)i in (2fy) at each iteration. 



4. The case of separable problems 

Let us now apply the method to the separable problem ([2]). By setting the (block) vectors 



Ik 



p = ( Pi > 



Pk 



q = { qi , ■■ 

one then obtains (see (fT2|)). 

q = e ® go + hA (8 I m p, p = e®po + hA® I m VJ7(q), 
i.e., since Ae — c= (ci, . . . , Ck) T , 

<\ = e®q + hc®p a + h 2 A 2 ®I m VU{(\). (25) 
Moreover, taking into account (|9|)- (fT2)) and (|17p . one obtains 

A 2 =l s _ 1 X s Vj„ i n. (26) 



The new approximations to q(to + h) and p(to + h) are then given by 

go + hp + h 2 e T flA ® I m VJ7(q), p + he T Q ® 7 m VE/(q), 

respectively. By using similar arguments as those given in Section [5] (see (HU), we set 

q = e® q a + hcgipQ + h 2 l s ^iX s ® I m 7, 

thus obtaining the following equation (which is the analogous of (|15[) ) : 

Firt) = 7 - (Vf-iSl ® ^m) Vf7 (e ® go + he ® p + h 2 I s -yX s ® I m j) = 0. (27) 

Similarly as what seen in Section [31 the application of the simplified Newton iteration for solving 
([27)1 then gives, by virtue of dTUJ) and ([T7]). and setting G = V 2 U(q ), 

[l-h 2 X 2 ®G Q ]A i = -F{ 1 e ), Y +1 =J e + & e , 1 = 0,1,..., (28) 

which, as in the previous case, has (block) size s, rather than k. The problem is then exactly that 
seen in (|18[) . via the formal substitutions 



h — > h 2 , X s — > X 2 . (29) 

This means that we can repeat similar steps for the blended solution of f28l) . by following the same 
arguments seen in Section [3[ In more details, ([T9)) -([23 |) can be repeated, by considering the formal 
substitutions ([2T))) and, moreover, 

P — > p 2 : hm — > Im- 
Also in this case 13, U|, the optimal choice of the parameter p turns out to be given by ([24]) . 



5. Numerical Tests 



We here consider a model problem to test the proposed algorithms and methods, in order to 
confirm the usefulness of the proposed approach. In particular, it is clear that a Newton-type 
iteration, like ([T8)) and (f28|) . works well when the linear part of the problem is significant. For this 
purpose, we consider the following polynomial Hamiltonian, 

H(q,p) = \v 2 - 10V - \q 2 -fcz+i), (30) 

from which we derive the following special second order problem, 

q" = 10 4 q (4q 3 - 3q 2 - 2q + l) , t G [0, 100], q(0) = 0, q'(0) = l. (31) 
For solving ([3"Tjl , we use the following fourth-order numerical methods: 

• the symplectic 2-stages Gauss method (GAUSS2); 

• HBVM(8,2) which is energy conserving, for the problem at hand. 
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Table 1: total number of iterations for solving the discrete problems with the specified stepsize h (- if no convergence). 



h 


secor 
blended 


GAl 

id order 
fixed-point 


JSS2 

firsl 
blended 


; order 
fixed-point 


secoi 
blended 


HBVr 

id order 
fixed-point 


,4(8,2) 

firsl 
blended 


; order 
fixed-point 


icr 3 
5 • i(r 3 

1(T 2 


664545 
242844 


690197 


952902 
308406 


1217673 


660317 
228242 
194163 


695765 
223883 


947618 
293949 
253049 


1225318 
424402 



For both methods, we consider a fixed-step implementation with stepsize h, with the generated 
discrete problems solved either with a fixed-point iteration or with a blended iteration, which have 
approximately the same cost, in terms of function evaluations. Moreover, we also compare the second 
order implementation described in Section^ with the equivalent first order Hamiltonian formulation 
of the problem, as described in Section [3] Table [1] summarizes the obtained results, in terms of total 
number of iterations (blended or fixed-point) for covering the specified integration interval. From 
the listed results, one deduces that the second order formulation of the problem is less demanding 
in terms of needed iterations. Moreover, the blended iteration turns out to be both more efficient 
and robust than the fixed-point iteration. 

Finally, in Figures [THD we plot the phase portraits of the numerical solutions, for the two methods 
and the various stepsizes, along with the corresponding error in the numerical Hamiltonian. As one 
can see, the phase portraits of the HBVM(8,2) method are always correct, whatever the used stepsize 
(see Figure Q] and the left plot in Figure [2]), since the Hamiltonian is conserved (up to round-off 
errors), as is shown in the right plot of Figure [5] This is not true for the GAUSS2 method, for which 
the error in the Hamiltonian depends on the used stepsize, as is shown in Figure [4j thus causing 
drawbacks in the corresponding phase portraits of the numerical solution, unless the stepsize is very 
small (see the two plots of Figure [3]) . 

From the numerical tests, one can then conclude that the proposed blended implementation of 
HBVMs turns out to be robust and efficient. Moreover, the energy-conserving property of such 
methods turns out to be very remarkable, with respect to standard symplectic methods. Finally, 
the second order formulation of HBVMs greatly improves their performance. 
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